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Abstract  ' 

A wavesplitting  procedure  is  proposed  by  which  noise  information  can  be  derived  from  CFD  calculations 
on  the  aeroengine  fan  stage.  Noise  propagation  in  the  ducted  regions  is  compared  with  well-understood 
linear  behaviour  in  parallel  wall  ducts.  Deviations  from  this  behaviour  are  used  to  highlight  important 
features  of  the  flow  solution.  These  include  genuine  flow  features  such  as  non-linear  acoustic  interaction 
as  well  as  dissipation  and  boundary  condition  errors  deriving  from  the  numerical  solution  of  the  equations. 
Properly  applied,  the  method  provides  quantitative  noise  source  amplitudes,  accounting  for  (modest) 
reflections  from  the  boundaries  of  the  CFD  domain.  At  the  same  time  confidence  can  be  gained  that  the 
CFD  results  are  accurate  for  the  wavelengths  and  frequencies  being  analysed,  and  that  the  CFD  domain 
sufficiently  covers  the  region  of  interest. 

Examples  are  given  of  how  the  method  can  be  applied  to  steady  and ‘'unsteady  CFD  calculations. 
Limitations  of  the  method  are  also  discussed.  ' 

Introduction 

Wavesplitting  (decomposition  of  a flow  field  into  upstream  and  downstream  travelling  eigenmodes)  is 
widely  used  in  CFD  boundary  conditions  to  determine  the  direction  in  which  information  is  travelling,  and 
hence  to  prevent  spurious  reflections  at  the  boundaries  of  the  domain.  Tyler  and  Sofrin  (1961)  first 
derived  the  acoustic  modes  for  uniform  axial  flow  in  cylindrical  and  annular  ducts.  This  formulation  has 
been  widely  and  successfully  applied  to  boundary  conditions:  Giles,  for  example,  used  the  modes  with 
varying  levels  of  approximation  to  produce  ID  and  2D  non-reflecting  conditions  (Giles  1990,  Saxer  and 
Giles  1993).  Other  authors  have  proposed  a complete  3D  formulation.  The  method  has,  however,  clear 
limitations.  Ducts  are  not  in  general  parallel  annulus.  Flow  is  generally  non-uniform  radially.  Amplitudes 
(particularly  with  regard  to  rotor  alone  tones)  are  so  high  that  the  acoustic  perturbations  behave  non- 
linearly  (Morfey  and  Fisher  1970).  Moreover,  in  addition  to  these  physical  features,  CFD  solutions  have 
their  own  characteristics:  smoothing,  grid  and  boundary  conditions  all  affect  the  noise  prediction  process. 
Rowley  and  Colonius  (1998)  developed  numerically  non-reflecting  boundary  conditions  by  calculating 
linear  modes  using  the  true  mean  flow  profile  and  included  in  their  analysis  the  discretisation  of  the 
particular  CFD  code  being  used.  - 

In  this  paper  a wavesplitting  method  is  developed  for  use  inside  the  computational  domain,  as  an  aid  to 
deriving  tone  noise  information  from  the  CFD  solution.  As  with  the  boundary  conditions,  a decision  has  to 
be  made  as  to  the  complexity  required.  ID  and  2D  analyses  are  insufficient  for  the  general  case.  Instead  a 
3D  uniform  axial  flow  method  based  on  Tyler  and  Sofrin  modes  is  adopted.  This  has  many  advantages 
over  more  complex  methods: 

• The  modes  can  be  calculated  quickly  using  standard  Bessel  functions. 

• The  modes  are  orthogonal,  making  decomposition  easier  and  faster. 

• The  behaviour  of  the  Tyler  and  Sofrin  modes  is  well  understood. 

• The  resulting  wavesplitting  routines  are  general  across  all  CFD  codes. 
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Two  methods  are  given  for  splitting  the  flow  field  at  a given  axial  plane  into  Tyler  and  Sofrin  modes, 
depending  on  whether  or  not  vortical  as  well  as  acoustic  waves  are  present  in  the  flow.  The  output  of  each 
method  is  a number  representing  amplitude  and  phase  for  the  upstream  and  downstream  travelling  acoustic 
waves  at  each  circumferential  and  radial  harmonic. 

Tyler  and  Sofrin  modes  propagate  either  at  constant  amplitude  (cut-on  waves)  or  with  log  amplitude 
varying  linearly  with  axial  distance  (cut-off  waves).  In  both  cases  the  phase  varies  linearly  with  distance. 
The  axial  wavenumber,  representing  the  rate  of  change  of  amplitude  and  phase,  can  be  calculated 
analytically  from  the  circumferential  wavenumber,  the  wave  frequency  and  the  (axial)  Mach  number  of  the 
mean  flow  (appendix  1).  Because  the  Tyler  and  Sofrin  modes  do  not  agree  exactly  with  either  the  physical 
or  the  CFD  calculated  modes  (for  the  reasons  outlined  above),  the  behaviour  of  the  calculated  wave 
amplitude  varies  from  this  ‘ideal’.  By  plotting  log  amplitude  and  phase  against  axial  distance  it  is  possible 
to  isolate  these  regions,  and  hence  to  direct  further  work  at  establishing  whether  the  cause  of  the 
discrepancy  is  physical  (eg  non-linearity,  hade  angle,  non-uniform  mean  flow)  or  numerical  (eg  dissipation 
error,  reflections  from  the  boundaries).  This  has  proved  an  essential  step  in  obtaining  reliable  noise 
information,  and  the  examples  below  show  how  the  analysis  might  proceed  from  the  basic  method  to 
derive  quantitative  estimates  of  source  noise. 

Example  1:  Steady  CFD,  Boundary  Conditions,  Non-linear  effects 

Before  describing  the  mechanics  of  the  method  it  is  useful  to  consider  an  example  to  illustrate  the 
progression  of  thought.  Fig  1 shows  the  result  of  a steady  CFD  calculation  of  a research  fan  blade  at  rig 
scale,  aimed  at  predicting  fan-forward  rotor  alone  noise.  The  pressure  has  been  Fourier  transformed 
circumferentially  to  decompose  it  into  the  harmonics  of  blade  passing  frequency  that  would  be  detected  by 
a stationary  microphone  at  the  outer  wall.  The  signal  amplitude  at  the  first  blade  passing  harmonic  is 
plotted  on  a log  scale  (in  fact,  as  Sound  Pressure  Level,  SPL1)  against  distance  upstream  of  the  fan  tip 
leading  edge.  The  large  oscillations  in  amplitude  are  known  from  measurements  to  be  false. 


0.8  0.6  0.4  0.2  0.0 

Axial  Distance  from  Fan  Tip  L/E  (Normalised  by  Fan  Diameter) 


Figure  1.  Steady  CFD  Solution  for  Fan  Rotor  Alone  Noise;  SPL  on  Outer  Wall 


1 Sound  Pressure  Level  (SPL)  is  defined  in  the  usual  way  as  201og)0(Prms/  2.10'5Pa) 
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The  starting  strength  close  to  the  fan  is  known  to  be  roughly  correct,  but  the  signal  (if  genuine)  is 
contaminated  by  oscillations  in  the  solution,  making  it  hard  to  quantify  the  noise  produced.  Two 
lengthscales  can  be  distinguished;  a short  wavelength  oscillation  with  about  eleven  peaks  over  the  axial 
domain,  and  a longer  wavelength  oscillation  with  troughs  at  x=0.5  and  x=0.3. 

The  first  step  in  the  analysis  procedure  is  to  split  the  pressure  signal  into  Bessel-Fourier  harmonics.  It  will 
be  seen  later  that  over  the  majority  of  the  domain  the  first  radial  harmonic  dominates  the  signal  at  the  outer 
wall.  This  harmonic  is  shown  as  the  long  dashes  in  fig  2,  compared  to  the  outer  wall  values  from  fig  1 . 
The  longer  lengthscale  oscillation  has  disappeared,  leaving  a smoothly  varying  version  of  the  shorter 
wavelength  disturbance.  Thus  the  longer  wavelength  disturbance  was  due  to  interference  at  the  outer  wall 
between  different  radial  harmonics. 


0.8  0.6  0.4  0.2  0.0 

Axial  Distance  from  Fan  Tip  L/E  (Normalised  by  Fan  Diameter) 


Figure  2.  Steady  CFD  for  Rotor  Alone  Noise,  Bessel-Fourier  Decomposition  and  Wavesplitting 


In  some  cases  this  level  of  analysis  is  sufficient.  The  radial  harmonic  breakdown  along  the.duct  gives  a 
great  deal  of  information  about  both  the  source  amplitude  and  the  accuracy  of  the  code,  as  will  be 
discussed  later.  In  this  case,  however,  a further  step  is  necessary  to  isolate  the  cause  of  the  short 
wavelength  oscillations.  The  periodic  nature  of  the  oscillations  suggests  some  form  of  standing  wave,  and 
this  particular  case  was  run  with  ID  non-reflecting  boundary  conditions  at  the  inlet  boundary,  which  are 
known  to  be  poor  for  the  wavelengths  and  frequencies  of  interest. 

By  decomposing  the  axial  velocity  also  into  Bessel-Fourier  harmonics  the  pressure  variation  represented 
by  the  long  dashes  in  fig  2 can  be  split  into  forward  and  backward  travelling  components.  The  short 
dashes  in  that  figure  represent  the  upstream  travelling  component.  The  result  is  a smooth  line, 
demonstrating  that  the  oscillations  observed  at  the  outer  wall  are  indeed  due  to  the  presence  of  a 
downstream  wave  reflected  from  the  inlet  boundary. 
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Axial  Distance  from  Fan  Tip  L/E  (Normalised  by  Fan  Diameter) 


Figure  3.  Steady  CFD  for  Rotor  Alone  Noise:  First  Three  Radial  Harmonics  of  Upstream  and 
Downstream  Waves 


Figure  3 shows  upstream  and  downstream  wave  components  of  the  noise  for  the  first  three  radial 
harmonics  at  the  first  blade  passing  frequency.  The  third  harmonic  is  cut-off  and  will  be  discussed  later. 
The  first  two  harmonics  are  cut-on,  and  Tyler  and  Sofrin’s  linear  method  would  predict  propagation  at 
constant  amplitude.  The  upstream  waves,  however,  decay  considerably,  especially  in  the  vicinity  of  the 
fan.  Fig  4 shows  the  first  radial  harmonic,  compared  with  Morfey  and  Fisher’s  ID  non-linear  decay  theory 
(Morfey  and  Fisher  1970).  The  agreement  is  not  expected  to  be  exact,  given  the  assumptions  implicit  in 
the  ID  calculation.  However,  the  relatively  close  agreement  near  the  fan  suggests  that  the  decay  is  in  the 
main  genuine  and  due  to  non-linear  effects.  The  less  good  agreement  further  from  the  fan  is  thought  to  be 
due  to  the  increasingly  coarse  grid  away  from  the  fan,  which  has  the  effect  of  attenuating  the  higher 
circumferential  harmonics. 
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Axial  Distance  from  Fan  Tip  L/E  (Normalised  by  Fan  Diameter) 


Figure  4.  Steady  CFD  for  Rotor  Alone  Noise,  Comparison  of  First  Radial  Harmonic  with  Morfey- 
Fisher  ID  Non-linear  Decay  Theory 

The  downstream  travelling  waves  (dashed  lines  in  fig  3)  are  lower  amplitude  and  are  less  affected  by  non- 
linear decay.  Hence  they  travel  at  more  nearly  constant  amplitude  over  the  majority  of  the  domain.  The 
high  ratio  of  reflected  to  incident  amplitude  at  the  inlet  to  the  domain  confirms  that  the  boundary 
conditions  are  indeed  behaving  poorly.  The  upstream  wave  amplitude,  however,  allows  a quantitative 
noise  estimate  to  be  obtained  without  contamination  from  these  reflected  waves,  provided  that  the  response 
of  the  fan  to  the  downstream  waves  can  be  ignored.  In  this  particular  case  the  non-linear  decay  is 
sufficient  to  mask  any  minor  changes  in  source  amplitude,  and  the  result  is  unlikely  to  be  affected  by  the 
presence  of  the  reflected  waves. 

The  oscillations  in  amplitude  of  the  reflected  waves  near  the  fan  leading  edge  are  a result  of  the 
assumptions  in  the  mathematics  of  the  wavesplitting.  The  relationship  between  pressure  and  velocity  as 
calculated  by  the  CFD  code  are  slightly  different  to  the  analytic  values,  due  to  discretisation  effects  and  the 
assumptions  inherent  in  the  analysis.  Hence  there  is  a limit  to  the  extent  to  which  upstream  and 
downstream  waves  can  be  discriminated.  This  can  be  a particular  problem  near  the  cut-off  boundary,  as 
explained  under'Bxample  2.  The  results  in  figure  3 show  that  upstream  and  downstream  waves  of  this  type 
can  be  discriminated  to  within  15-20dB  of  the  higher  signal  amplitude.  ' 

The  third  harmonic  shown  in  figure  3 is  cut-off,  and  decays  rapidly  with  axial  distance  to  a ‘floor’  level  of 
around  llOdB.  This  floor  level  is  set  by  ‘leakage’  between  radial  harmonics.  The  actual  propagating 
modes  in  the  calculation  are  not  pure  Bessel-Fourier  harmonics,  both  because  of  discretisation  effects  and 
because  the  mean  flow  is  not  uniform.  Hence  the  Bessel-Fourier  breakdown  in  the  wavesplitting  process 
calculates  some  content  at  all  radial  harmonics  even  if  only  one  mode  is  present.  Nonetheless,  the 
discrimination  between  radial  modes  (over  40dB  from  the  first  harmonic  amplitude  to  the  ‘floor’  level  in 
fig  3)  is  more  than  sufficient  for  most  purposes. 

The  reason  for  the  poor  discrimination  between  upstream  and  downstream  waves  near  the  domain  exit  at 
this  third  harmonic  is  not  known.  Typically  (away  from  the  cut-off  point  itself)  reasonable  discrimination 
is  obtained  for  both  cut-on  and  cut-off  waves. 

Figure  5 shows  the  amplitude  of  the  upstream  component  of  the  third  harmonic  compares  well  with  the 
exponential  decay  predicted  in  Tyler  and  Sofrin’s  theory. 
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Axial  Distance  from  Fan  Tip  L/E  (Normalised  by  Fan  Diameter) 

Figure  5.  Steady  CFD  for  Rotor  Alone  Noise,  Comparison  of  Third  Radial  Harmonic  Upstream 
Travelling  Wave  with  Linear  Theory 


The  noise  output  of  the  fan  can  only  be  estimated  in  the  region  where  the  mesh  is  sufficiently  fine.  In  this 
case,  the  noise  output  (fig  3)  is  163dB  SPL  in  the  first  radial  harmonic,  152dB  in  the  second,  and 
negligible  output  in  the  third  and  higher  harmonics,  all  at  a position  0.2  diameters  upstream  of  the  fan. 
Note  that  this  figure  could  not  easily  have  been  obtained  from  either  the  wall  pressure  or  the  harmonic 
breakdown  of  pressure  alone  (fig  2). 

Figure  4 suggests  that  there  is  still  significant  non-linear  behaviour  upstream  of  the  0.2m  point.  Depending 
on  the  use  to  which  the  data  is  to  be  put,  it  may  be  important  to  contain  the  non-linear  behaviour  within  the 
CFD  domain.  In  this  case  a second  CFD  calculation  would  be  required  with  increased  mesh  density  in  the 
upstream  region. 

To  summarise,  in  this  first  example  the  wavesplitting  method  has  been  used 

• to  demonstrate  that  the  observed  oscillations  in  pressure  are  due  to  reflections  from  the  inlet  boundary 
and  interference  between  radial  harmonics, 

• to  match  the  behaviour  of  the  CFD  solution  to  expectation  and  hence  to  give  confidence  in  the 

solution,  " 

"•  to  gauge  the  significance  of  non-linear  effects,  and 

• to  obtain  a quantitative  measure  of  the  fan-forward  rotor  alone  noise  in  the  presence  of  significant 
reflections  from  the  upstream  plane. 

Method 

The  post-processing  required  for  the  wavesplitting  method  can  be  broken  down  into  the  following  steps: 

1 . Fourier  transform  in  time  (non-linear  unsteady  CFD  codes  only). 

2.  Interpolate  onto  axial  planes. 

3.  Decompose  into  Bessel-Fourier  harmonics. 

4.  Calculate  linear  propagation  parameters  at  each  axial  plane  based  on  average  mean  flow. 

5.  Use  linear  theory  to  split  Bessel-Fourier  harmonics  into  upstream  and  downstream  components. 

It  can  be  useful,  in  comparing  the  relative  strengths  of  the  radial  harmonics,  to  add  a sixth  step: 

6.  Calculate  acoustic  energy  flux  for  upstream  and  downstream  components. 
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In  linear  theory  with  uniform  axial  mean  flow  and  parallel  annulus  ducts  (Tyler  and  Sofrin  1961),  cut-on 
modes  propagate  with  uniform  amplitude  and  cut-off  modes  with  exponential  decay.  Both  types  propagate 
with  phase  varying  linearly  with  axial  distance.  The  decay  rates  and  axial  wavenumber  are  calculated  at 
step  4,  and  so  by  plotting  amplitude  on  a log  scale  (eg  as  SPL)  and  phase  on  a linear  scale,  it  is  easy  to 
compare  the  modal  behaviour  with  theory.  Variations  from  linear  behaviour  can  then  be  easily  identified 
and  subjected  to  further  analysis,  as  demonstrated  in  the  examples  above  and  below. 

Step  1.  Fourier  transform  in  time. 

Both  the  theoretical  linear  calculations  and  the  wavesplitting  require  data  at  fixed  position  and  fixed 
frequency.  In  steady  CFD  calculations  on  rotor  blades  the  frequency  is  fixed  by  the  circumferential 
harmonic  (that  is,  a first  blade  passing  harmonic  in  the  circumferential  direction  will  be  heard  at  first  blade 
passing  frequency  by  a stationary  observer).  In  linear  unsteady  calculations  the  frequency  is  fixed  in  the 
frame  of  reference  of  the  CFD  domain,  and  again  the  frequency  in  the  stationary  domain  can  be  easily 
calculated  for  a given  circumferential  harmonic.  In  non-linear  calculations  a Fourier  transform  in  time  is 
required  first  to  obtain  data  at  a given  frequency.  Such  a transform  is  easily  calculated  given  a periodic 
converged  solution. 

Step  2.  Interpolate  onto  axial  planes. 

Some  care  is  required  in  this  step,  particularly  in  the  presence  of  shocks,  where  a simple  linear 
interpolation  can  lead  to  shock  smoothing.  In  cases  with  structured  or  semi-structured  grids,  greater 
accuracy  can  sometimes  be  obtained  by  applying  a circumferential  FFT  first  (that  is,  to  bring  the  first  part 
of  step  3 forward),  and  then  to  interpolate  axially  for  amplitude  and  phase  at  a given  plane. 

Step  3.  Decompose  into  Bessel-Fourier  harmonics. 

This  process  consists  of  a Fourier  transform  circumferentially,  followed  by  a Bessel  function 
decomposition  in  the  radial  direction.  The  latter  is  aided  by  the  orthogonality  condition  (appendix  1).  The 
choice  of  flow  variables  to  be  decomposed  depends  on  the  type  of  wavesplitting  technique  to  be  employed, 
as  described  under  step  5. 

Step  4.  Calculate  linear  propagation  parameters  at  each  axial  plane  based  on  average  mean  flow. 

The  linear  propagation  parameters  (axial  wavenumbers  for  upstream  and  downstream  propagating  waves) 
are  calculated  for  each  circumferential  and  radial  harmonic,  based  on  a duct  with  parallel  annulus  walls 
and  uniform  axial  mean  flow  (appendix  1).  Given  the  approximations  of  the  method  a simply  calculated 
Mach  number  based  on  area  average  axial  velocity  and  temperature  (or  equivalent)  at  each  axial  plane  is 
sufficient  to  define  the  equivalent  mean  flow. 

Step  5.  Use  linear  theory  to  split  Bessel-Fourier  harmonics  into  upstream  and  downstream 
components. 

This  step  is  not  always  necessary.  In  particular  if.  the  boundary  conditions  in  the  CFD  calculation  are 
known  to  be  genuinely  non-reflective,  and  the  duct  in  which  the  analysis  is  being  performed  is  fairly 
uniform,  then  the  noise  can  be  assumed  to  be  travelling  in  a single  direction  only,  and  the  pressure 
information  alone  is  sufficient.  In  many  cases,  however,  the  effect  of  the  boundary  conditions  and  duct 
geometry  is  not  known  at  acoustic  wavelengths  and  frequencies,  and  this  step  is  required. 

Two  methods  of  performing  the  wavesplit  are  detailed  in  appendix  1.  The  first,  simpler  method  uses  just 
the  pressure  and  the  axial  velocity  from  the  CFD  solution,  and  is  based  on  the  assumption  of  irrotational 
flow  (that  is  in  total,  not  just  irrotational  mean  flow).  This  is  often  a good  approximation  in  the  inlet  duct. 
It  is  unlikely  to  be  sufficient  in  the  OGV  exit  duct,  or  in  other  calculations  where  vortical  waves  are 
prominent  in  the  flow  solution.  Boundary  conditions  too  (even  if  termed  non-reflective)  can  introduce 
vorticity  at  the  inlet  plane. 

The  second  method  is  still  based  on  the  assumption  of  uniform  axial  mean  flow,  but  allows  for  vortical 
waves  in  the  solution.  This  method  is  more  complex,  using  the  divergence  of  the  flowfield  to  remove  the 
effects  of  vorticity. 

Step  6.  Calculate  acoustic  energy  flux  for  upstream  and  downstream  components. 

Usually  the  most  useful  quantitative  output  from  the  analysis  method  is  the  noise  amplitude  (SPL)  in  the 
different  circumferential  and  radial  harmonics.  Sometimes,  however,  it  is  useful  to  compare  the  noise 
generated  in  the  different  harmonics,  or  to  combine  the  results  into  a single  figure.  Tester  (1972)  gives  a 
formula  (reproduced  in  appendix  1)  for  acoustic  energy  flux  based  on  the  same  uniform  axial  flow  and 
parallel  duct  assumptions  as  the  theoretical  calculations  in  step  4 and  the  irrotationality  assumption  of  the 
first  method  in  step  5.  The  calculation  can  be  applied  directly  to  the  interpolated  flow  variables  at  each 
axial  plane  at  step  2,  but  it  can  also  be  applied  to  the  upstream  and  downstream  components  of  the 
individual  circumferential  and  radial  harmonics.  Under  the  given  assumptions  these  behave  independently 
for  cut-on  modes  (that  is,  the  total  energy  flux  is  the  sum  of  the  constituents).  For  cut-off  modes  energy  is 
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transported  by  combinations  of  upstream  and  downstream  propagating  waves,  and  thus  the  calculation  has 
to  be  performed  prior  to  the  wavesplit. 

By  summing  the  energy  flux  over  all  the  modes  in  one  direction  only  (upstream  or  downstream),  it  is 
possible  to  calculate  an  overall  source  level  in  the  cut-on  modes  disregarding  any  waves  reflected  from  the 
boundary  of  the  CFD  domain. 

Example  2:  Unsteady  CFD,  Boundary  Conditions,  Dissipation  Errors 

This  example  comes  from  an  investigation  into  non-reflective  boundary  conditions.  Two  separate  3D 
unsteady  linear  calculations  were  performed  in  a parallel  annulus  duct  with  two  types  of  boundary 
condition  at  the  inlet  (Giles  type  ID  and  2D  non-reflecting,  Giles  1990,  Saxer  and  Giles  1993).  A uniform 
rectangular  grid  was  used,  with  61  points  axially  and  41  points  circumferentially  and  radially, 
representative  of  a short  section  of  the  inlet  duct  for  a typical  aeroengine  fan  at  rig-scale  (Outer  diameter 
around  0.43m).  The  base  flow  for  the  linear  calculation  was  set  to  be  uniform  axial  flow,  again 
representative  of  a fan  inlet  duct  at  transonic  conditions.  In  each  case  a 26-lobed  upstream  travelling 
acoustic  wave  was  input  at  the  outflow  boundary,  and  the  output  analysed  using  the  method  described  in 
the  previous  section. 

Figure  6 shows  some  of  the  results  from  the  calculation  with  ID  non-reflecting  boundary  conditions.  The 
input  wave  had  constant  amplitude  radially.  This  implies  content  at  all  radial  harmonics,  and  so  it  is 
unsurprising  that  the  wall  amplitude  (solid  line  in  Figure  6)  varies  considerably  along  the  length  of  the 
duct.  The  total  amplitude  of  the  first  radial  harmonic  (dashed  line)  also  varies,  due  to  reflections  from  the 
inlet  boundary.  When  the  wavesplittirig  process  is  applied  using  the  static  pressure  and  axial  velocity 
(dotted  line),  the  long  lengthscale  oscillations  are  removed,  but  more  oscillations  are  present  at  a shorter 
lengthscale. 

The  oscillations  in  wave  amplitude  are  the  result  of  vorticity  introduced  by  the  boundary  conditions  at  the 
inlet  boundary.  The  reflected  vortical  waves  interfere  with  the  axial  velocity  field,  preventing  an  accurate 
analysis  of  the  acoustic  waves.  This  effect  can  be  removed  by  applying  the  wavesplit  method  using  static 
pressure  and  the  divergence  of  the  velocity  field  (see  appendix  1).  The  result  of  this  analysis  is  shown  as 
the  dash-dot  line  in  Figure  6.  All  oscillations  have  been  removed,  leaving  a uniformly  decaying  amplitude 
from  which  a quantitative  estimate  can  be  made  of  the  dissipation  error  inherent  in  the  numerical 
calculation  ( 1 5dB  per  duct  diameter).  Similar  calculations  could  be  performed  using  the  phase  of 
individual  waves  to  determine  the  propagation  phase  error. 


Axial  Distance  (Normalised  by  Duct  Diameter) 


Figure  6.  Unsteady  CFD,  Wavesplit  Based  on  Ua  and  div.U 
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Figure  7.  Unsteady  CFD:  Upstream  and  Downstream  Waves 

Full  results  from  the  wavesplitting  procedure  at  the  first  and  second  radial  harmonics  are  presented  in 
Figure  7.  The  second  harmonic  is  very  close  to  cut-off  (cut-off  ratio  1.01),  and  illustrates  some  of  the 
difficulties  commonly  encountered  close  to  the  cut-on/cut-off  boundary.  The  upstream  wave  shows  an 
exceptionally  fast  rate  of  decay  with  axial  distance,  and  the  downstream  wave  in  particular  shows 
significant  non-linear  behaviour,  in  that  the  (log)  amplitude  varies  non-linearly  with  axial  distance. 

Because  information  in  near  cut-off  modes  travels  slowly  along  the  duct  there  is  more  time  for  genuine 
non-linear  and  viscous  phenomena  to  affect  noise  propagation  between  axial  planes.  The  present 
calculations,  however,  were  performed  with  an  inviscid  linear  code,  and  so  neither  of  these  factors  can 
explain  the  high  decay  rate  of  the  upstream  wave.  The  decay  here  is  a result  of  the  numerical  smoothing 
applied  in  the  CFD  code. 

In  a CFD  calculation  smoothing  is  balanced  by  the  rate  at  which  the  wave  information  travels  upstream. 
Hence  waves  near  cut-off  tend  to  be  heavily  attenuated.  The  principle  can  illustrated  using  a simplistic  2D 
wave  equation  solved  analytically  with  second  order  smoothing  on  the  pressure  term.  Using  the 
nomenclature  defined  in  the  appendix, 

du  *7 

_ = -cvp 

dt 

— = -cVai  + eV3/? 

dt 

The  level  of  smoothing  £ is  assumed  to  be  small  and  fixed  independently  of  the  frequencies  and 
wavelengths  present  in  the  solution.  For  a cut-on  upstream  plane  wave  at  a given  frequency  CO  it  is 
possible  to  show  (with  some  algebra)  that  to  first  order  in  £the  axial  wavenumber  of  the  smoothed  wave  is 
X = X +iem2/2c2uinf 

where  winf  = - c2X/(0  is  the  rate  at  which  information  travels  upstream,  and  X is  the  axial  wavenumber  in 
the  absence  of  smoothing.  Thus  the  rate  of  decay  with  axial  distance  (determined  by  the  imaginary  part  of 
the  axial  wavenumber)  is  £C02/2c'Mmr  . Hence  waves  near  the  cut-off  boundary  (wjnf  — >0)  are  heavily 
attenuated. 

In  addition  to  the  greater  non-linear  and  smoothing  effects,  the  wavesplitting  process  itself  is  liable  to 
errors  close  to  the  cut-on/cut-off  boundary.  At  this  boundary  the  upstream  and  downstream  linear  waves 
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coalesce.  In  theory  this  could  give  large  numerical  errors  in  the  wavesplitting  solution  (mathematically, 
fm  and  fds  get  very  close  in  equations  18,  19  of  the  appendix)2.  In  practice  this  is  rarely  a problem:  even 
with  a cut-off  ratio  of  1.01  there  is  sufficient  discrimination  to  prevent  large  numerical  errors. 

Of  much  greater  significance  is  the  fact  that  the  discrepancies  between  the  CFD  calculated  modes  and  the 
analytical  modes  grow  large  near  the  cut-on/cut-off  boundary.  In  the  analytic  axial  mean  flow  formulation 
the  axial  wavenumber  X becomes  very  sensitive  to  flow  conditions,  and  in  particular  to  the  cut-off  ratio 
S,  =o)/c)i(l  -U1)'  Indeed,  it  can  be  seen  from  equation  7 of  the  appendix  that  at  the  boundary  itself  the 

derivative  of  X with  respect  to  (;  becomes  infinite.  It  is  unsurprising  therefore  that  the  analytic  modes  and 
those  calculated  by  CFD  (with  non-uniform  mean  flow,  smoothing  and  discretisation  errors)  can  be  very 
different  near  the  boundary. 

The  following  graph  illustrates  the  problems  that  can  occur  near  the  cut-off  boundary.  The  simple  case  is 
taken  of  a 2D  wave  equation  solved  using  analytic  derivatives  in  time  and  axially  (x  direction),  but  a 
simple  central  difference  numerical  derivative  in  the  y direction.  Forty  points  per  wavelength  were  used, 
which  is  usually  sufficient  to  give  reasonable  wave  propagation.  An  exact  modal  solution  to  the  numerical 
problem  was  calculated  for  a range  of  frequencies,  to  which  the  wavesplitting  technique  was  applied.  The 
graph  shows  the  amplitude  of  the  forward  and  reverse  waves  calculated  by  the  method  for  points  close  to 
the  cut-off  boundary.  As  expected,  the  errors  grow  large  close  to  the  boundary,  and  very  close  to  the 
boundary  the  calculated  waves  in  both  directions  can  be  much  larger  in  amplitude  than  the  genuine  wave. 
Note,  though,  that  the  errors  are  highly  concentrated  in  the  region  of  the  boundary:  away  from  this  point 
(towards  either  cut-on  or  cut-off)  the  errors  are  small  and  the  wavesplitting  technique  is  successful. 


Cut— on  Ratio 

Figure  8.  Errors  in  Wavesplitting  Method  Close  to  the  Cut-Off  Boundary 


2 With  the  second  wavesplitting  method,  based  on  pressure  and  the  divergence  of  the  velocity,  numerical 
errors  could  also  cause  a problem  in  regions  of  low  or  zero  mean  flow.  This  is  because,  although  the  axial 
wavenumbers  of  the  upstream  and  downstream  waves  remain  discrete  with  decreasing  mean  flow,  the 
wavesplitting  parameters  /“'  and  fds  (equation  22  of  the  appendix)  become  identical. 
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Returning  to  the  CFD  example,  Figure  9 shows  the  results  of  the  wavesplitting  analysis  for  the  ID  and  2D 
non-reflective  boundary  conditions.  The  ID  boundary  conditions  give  a reflected  wave  with  amplitude 
1 1 ,5dB  lower  than  the  incident.  The  second  radial  harmonic,  however,  shows  a reflection  only  3dB  lower 
than  the  incident  (accepting  the  errors  discussed  above  in  the  wavesplitting  for  this  harmonic).  Care  has  to 
be  applied  here,  in  that  boundary  conditions  can  scatter  an  incident  wave  at  one  radial  harmonic  into 
reflected  waves  at  other  harmonics.  However,  a more  detailed  analysis  of  the  ID  non-reflective  boundary 
conditions  used  in  this  case  has  shown  that  the  downstream  wave  at  the  second  radial  harmonic  is  indeed  a 
reflection  of  the  same  harmonic,  rather  than  a scattered  reflection  from  the  first  harmonic.  The  2D  non- 
reflective  boundary  conditions  give  a reflection  of  the  first  harmonic  30dB  lower  than  the  incident.  This 
level  is  close  to  the  level  at  which  forward  and  backward  waves  can  be  discriminated,  hence  the  oscillation 
with  axial  distance.  The  second  harmonic  is  reflected  at  a level  4dB  lower  than  the  incident.  A more 
detailed  analysis  of  the  boundary  conditions  has  shown  that  this  value  would  be  7.6dB  without  scattering 
from  the  higher  amplitude  first  radial  harmonic. 

In  both  cases  the  high  reflection  coefficients  for  the  second  harmonic  mirror  the  difficulty  in  wavesplitting 
close  to  the  cut-off  boundary. 
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Figure  9.  Unsteady  CFD:  Comparison  of  ID  and  2D  Boundary  Conditions 


Further  Example 

The  wavesplitting  method  has  been  applied  to  a flat  plate  wake/vane  interaction  case  (Wilson  2001).  In 
that  work  the  method  was  also  extended  to  provide  fully  3D  non-reflecting  boundary  conditions  at  inflow 
and  exit. 

Limitations 

The  major  limitations  of  the  analysis  procedure  have  been  identified  as  follows: 

1)  The  method  as  it  stands  is  unsuitable  for  use  in  regions  where  the  propagating  modes  are  known  to 
differ  strongly  from  Tyler  and  Sofrin  modes.  Examples  are: 

a)  Fan  exit,  where  the  high  degree  of  swirl  significantly  affects  the  propagating  modes. 

b)  Lined  ducts  (although  the  method  can  still  be  applied  in  any  hardwall  regions  upstream  and 
downstream  of  the  liner). 

2)  The  method  gives  poor  results  close  to  the  cut-on/cut-off  boundary  (as  discussed  under  example  2). 
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3)  The  wavesplitting  method  based  on  pressure  and  axial  velocity  is  limited  to  regions  of  irrotational 
flow.  This  can  be  a good  approximation  in  the  fan  inlet  region,  subject  to  the  boundary  conditions  at 
the  inlet  plane  producing  negligible  vorticity.  This  limitation  can  be  overcome  by  using  the  second 
wavesplitting  method,  based  on  pressure  and  the  divergence  of  the  velocity. 

Summary 

A procedure  has  been  developed  for  deriving  tone  noise  information  from  CFD  solutions  in  the  fan  stage. 
A wavesplitting  method  is  used  at  a number  of  axial  locations  upstream  of  the  fan  to  derive  upstream  and 
downstream  wave  information  from  the  CFD  calculated  flowfield.  The  development  of  these  waves  with 
axial  distance  is  then  compared  to  well-understood  linear  behaviour. 

The  method  is  equally  applicable  to  the  bypass  duct  downstream  of  the  outlet  guide  vane  (where,  like  the 
inlet,  the  mean  flow  is  approximately  axial).  Application  to  other  areas  is  subject  to  the  limitations 
outlined  in  the  previous  section. 

Examples  have  been  given  of  how  the  method  can  be  used 

• to  identify  regions  of  non-linear  behaviour 

• to  identify  and  eliminate  the  effect  of  modest  spurious  reflections  at  the  domain  boundaries 

• to  identify  and  eliminate  the  effect  of  modest  vortical  waves  present  in  the  solution 

• to  quantify  numerical  dissipation  errors 
and  hence 

• to  quantify  tone  noise  generation. 
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Appendix  1 - Mathematics  of  Wavesplitting  Technique 


Nomenclature 

U Mean  flow  speed  (assumed  axial)  normalised  by  speed  of  sound 
p Mean  flow  pressure 

P Mean  flow  density 

c Speed  of  sound  based  on  mean  flow  pressure  and  density 
u Perturbation  velocity  normalised  by  speed  of  sound 
p Perturbation  pressure,  normalised  by  yp 

y Ratio  of  specific  heats  (assumed  constant) 

Wave  Solutions  of  the  Euler  Equations  in  Parallel  Annulus  Ducts 


The  linearised  Euler  equations  for  a perfect  gas  with  constant  ratio  of  specific  heats  and  uniform  mean 
flow  in  the  x direction  can  be  expressed  as 


1 9u  r,3u  _ 
-—  + U—-  = -Vp 
c dt  dx 


( 1) 


i^P 


-^  + U^  = -V.  u 

c dt  dx 


(2) 


Together  with  the  equation  of  state  and  the  energy  equation  this  represents  five  equations  in  five 
unknowns.  The  wave  solutions  (that  is,  solutions  for  the  flow  variables  that  preserve  their  shape  axially) 
are  well  known,  and  correspond  to  two  acoustic  waves,  two  vorticity  waves  and  a single  entropy  wave. 


For  a parallel  wall  annular  duct  the  acoustic  waves  can  be  written  as 


u = V<t.  ( 3 ) 

<*>  = ip/fa/c -Uk)  (4) 

where  p = pm  t w exp/(cor  - he-  mQ  )B‘  (r)  ,5, 


and  5*  (r)  is  a combination  of  the  J-  and  Y-type  Bessel  functions: 


KirhoJ^rf+bY.fcr),  (6) 

3(t) 

where  a,  b and  p*  are  chosen  to  satisfy  the  boundary  condition  ur  = — — = 0 on  both  the  inner  and  outer 

dr 

walls.  If-no  routine  is  available  to  give  these  functions,  a method  is  outlined  later  which  requires  only 
evaluation  of  the  individual  Bessel  functions. 

For  each  circumferential  harmonic  m and  radial  harmonic  k there  are  two  possible  axial  wavenumbers  A 
defined  by 


? -Uu/c±^ (oi/c)2  - n**(l -U2) 

1 -u1 

X is  either  real  or  complex  depending  on  the  sign  of  the  term  under  the  square  root 
consider  the  two  cases  separately: 

Cut-on  waves  (A.  real) 

It  can  be  seen  from  equation  5 that  if  X is  real,  then  the  wave  propagates  axially  at  constant  amplitude  and 
varying  phase.  For  forward  subsonic  mean  flow  0 < U<  1,  the  negative  root  in  equation  7 corresponds  to 
an  upstream  travelling  wave,  whilst  the  positive  root  corresponds  to  a downstream  wave.  Note  that  with 
mean  flow  U > 0 it  is  possible  for  X to  remain  negative  even  for  the  downstream  travelling  wave.  In  this 


(7) 

. It  is  convenient  to 
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case  the  wavefront  is  angled  in  the  same  direction  as  the  upstream  wave  and  indeed  appears  to  move 
upstream.  The  group  velocity,  however,  defining  the  direction  in  which  information  is  transported,  is  still 
downstream,  as  is  the  transport  of  energy. 

Cut-off  waves  (A  complex) 

If  A is  complex,  then  the  wave  decays  exponentially  in  the  direction  of  propagation.  For  forward  subsonic 
mean  flow  0 < U < 1,  the  root  of  equation  7 with  positive  imaginary  part  represents  an  upstream  wave 
decaying  with  upstream  distance,  and  the  root  with  negative  imaginary  part  represents  a downstream  wave 
decaying  with  downstream  distance.  The  phase  variation  is  fixed  by  the  first  term  in  equation  7,  which  is 
independent  of  m and  k.  Hence  all  cut-off  waves  have  the  same  ‘spiral  angle’  (that  is,  the  same  rate  of 
change  of  phase  with  axial  distance). 

In  each  case,  the  two  roots  of  equation  7 correspond  to  one  upstream  and  one  downstream  wave.  Entropy 
and  vorticity  waves  make  no  contribution  to  the  pressure,  and  so  at  a given  frequency  co,  circumferential 
mode  m and  radial  harmonic  k,  the  entire  unsteady  pressure  field  can  be  written  as  the  sum  of  the 
contributions  from  the  two  acoustic  waves: 

p = l pu'  (*)+  p'1'  (ar)]exp/((0?  - m9  )S*  (r),  where  ( 8 ) 


p“'  =c.ft.exp(-/A"'jc) 
p,h  = c.rf.exp(-  /A‘,'.v) 


( 9,  10 ) 


Thus  if  pm  and pds  are  known  at  any  single  axial  plane,  then  their  (theoretical)  values  can  be  calculated  at 
any  other  plane.  Two  methods  are  given  below  for  splitting  the  unsteady  pressure  into  its  upstream  and 
downstream  components,  depending  on  whether  or  not  the  flow  can  be  assumed  irrotational. 


Wavesplitting  in  Irrotational  Flow 

The  velocity  and  pressure  perturbations  of  individual  acoustic  waves  are  linked  by  equations  3 and  4. 
Entropy  waves,  furthermore,  make  no  contribution  to  either  pressure  or  velocity.  Hence  for  each  harmonic 
(CO,  m,  k)  the  axial  velocity  can  be  derived  from  the  upstream  and  downstream  components  of  pressure  (equ 
8)  as  follows: 

m,  = [»"'  (x)+  i<f  (x)Jexp/(cot  - w9  )Bkm  (r) , where  (11) 


u“'  = 

= /UV' 

/"'  = \u'/(co/c  - UT' ) 
/"'  =A‘,7(co/c-fA‘1') 


( 12,  13,  14,  15) 


Of  more  interest  is  the  inverse  process,  whereby  the  upstream  and  downstream  components  can  be 
determined  from  the  static  pressure  and  axial  velocity.  If  at  a given  axial  plane  and  these  two  variables  are 
decomposed  into  Bessel-Fourier  harmonics,  with  coefficients p\u\  then 


P = Pm  + P 

ux=tC+u?=rP“'+fd'P* 


Hence 

p-=(f-p 

p"! =(r  P -u,)Kr -r) 


( 18, 19) 
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Wavesplitting  in  the  Presence  of  Vortical  Waves 

In  the  presence  of  vortical  waves  the  previous  wavesplitting  method  fails,  because  the  axial  velocity  is  no 
longer  a function  of  the  acoustic  waves  alone.  Hence  another  method  is  presented  here,  based  on  pressure 
and  the  divergence  of  the  velocity  field  ( V • u ).  Note  that  this  method  too  is  based  on  the  assumption  of  a 
uniform  axial  mean  flow.  Thus,  whilst  small  amounts  of  vorticity  can  be  handled  using  this  method,  it  is 
not  suitable  for  highly  vortical  flows  such  as  those  found  between  an  aeroengine  fan  rotor  and  outlet  guide 
vane. 

The  divergence  of  the  velocity  for  acoustic  waves  can  be  calculated  from  the  two  pressure  coefficients 
(equs  9,  10)  using  equ  4; 

V-u‘'=VV 

=-(r!+My)r 
=-(r!+n:iV-/Wc-w) 

where  f^."‘  + (l*  ffco/c  - U\“  ) (21) 

With  some  manipulation,  this  can  be  reduced  to 

/"'  =-f(w/c-t/X“'),  (22) 

and  similarly  for  the  downstream  waves. 


(20) 


Both  V - u and  p are  identically  zero  for  vorticity  and  entropy  waves.  Hence  the  upstream  and  downstream 
components  can  be  obtained  from  V • u and  p in  an  identical  manner  to  the  irrotational  case  (equs  16  to 
19),  but  with  V - u instead  of  and  the  new  definitions  of  fm , fds  . 


V - u can  be  calculated  directly  from  the  CFD  solution  using  a simple  difference  sum.  A refinement  is  to 
use  the  flow  equations  (1, 2)  to  remove  the  axial  derivative  of  the  velocity,  such  that  the  wavesplit  can  be 
performed  using  only  data  from  a single  axial  plane.  The  following  analysis  is  for  isentropic  flow, 
although  the  result  also  holds  in  the  presence  of  entropy  waves. 

For  isentropic  flow,  given  the  normalisation  used,  the  continuity  equation  (2)  can  be  expressed  in  terms  of 
pressure: 


1 dp 
c dt 


T dp  „ 15  dur 

+ U—  = -V.u  = - V .u 


dx 


dx 


(23) 


2d  i 3(rwr ) 

where  V .u  = 

r dr 


, 1 due 
r dd 


(24) 


But  from  the  axial  component  of  the  momentum  equation  (1), 

dp=  i <K  u<K 

dx  c dt  dx 

Substituting  this  into  equ  23,  and  applying  a transform  in  the  time  direction  — — > iw, 

dt 


(25) 
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iw  dux 

— u*~u^r 

c ox 
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V dl,< 

- V .u *- 

dx 


^=-^.u+^(P-uux)y(i-u') 

V.u  = V .u  + ^ = -^C/2  V .u  + — (p  - Uux  )jy/(l  - U- ) 


(26) 


Thus  to  calculate  V • u from  data  at  a single  axial  location,  it  is  necessary  only  to  calculate  the  radial  and 
circumferential  derivatives  in  equ  24  and  to  substitute  the  result  together  with  values  for p and  w,  into  the 
equation  above.  Note  that  the  circumferential  derivative  is  trivial,  as  the  data  has  previously  been  Fourier 
transformed.  In  the  radial  direction  a simple  finite  difference  method  is  often  sufficient. 

The  use  of  a simple  finite  difference  sum  and  the  inclusion  of  the  analytic  flow  equations  do  make  the 
method  somewhat  approximate.  Nonetheless,  provided  the  vortical  waves  are  of  similar  amplitude  to  the 
acoustic  waves  (or  smaller)  it  is  effective  in  isolating  the  acoustic  component,  as  illustrated  by  example  2 
of  the  main  paper. 


Acoustic  Energy  Flux 


Tester  (1972)  gives  an  equation  for  modal  energy  flux  with  uniform  axial  mean  flow  (based  on  that  of 
Cantrell  and  Hart,  1964).  In  the  current  nomenclature,  at  any  given  frequency. 


( xmk 


p c-B 


mk 


Re  + Mxpmkp’nk  + M]pmkumk  +u,nku 


(27) 


where  pmk  and  umk  represent  coefficients  of  pressure  and  axial  velocity  at  a single  circumferential  (m)  and 
radial  ( k ) harmonic  at  the  given  frequency.  The  term 

Blk=\Bkm{r)\MT  (28)  - 

has  been  added  to  integrate  the  energy  flux  over  the  area  of  the  duct  (note  that  B7mk  = j B^(r)2 .2jirdr  is 
purely  real  in  this  equation). 

If  the  (m,k)  mode  is  cut-on,  then  the  upstream  and  downstream  travelling  waves  are  orthogonal  with 
respect  to  the  above  equation.  That  is,  the  total  energy  flux  in  the  mode  is  the  sum  of  that  calculated  from 
just  the  upstream  and  downstream  travelling  components.  If  the  mode  is  cut-off,  then  the  upstream  and 
downstream  travelling  waves  interact  to  give  energy  flux,  and  the  calculation  has  to  be  performed  prior  to 
the  wavesplit. 

Given  the  orthogonality  of  the  radial  harmonics,  the  total  acoustic  energy  flux  in  a given  circumferential 
harmonic  m can  also  be  calculated  by  integrating  the  radial  profiles  directly: 

J Rej/vC  +M,Pm Pi  +M)pmum  +uniu'mM ,}2itrdr  ( 29) 


where  pm  and  um  now  represent  radial  profiles  of  pressure  and  axial  velocity  components  at  the  m,h 
circumferential  harmonic.  Provided  that  a single  value  of  Mx  is  used  (rather  than  a radial  profile)  the  result 
is  equal  to  the  sum  of  equ  27  over  all  radial  harmonics  k. 


Calculation  of  Bkm(r) 

If  no  routine  is  available  to  give  these  functions,  a method  is  outlined  here  which  requires  only  values  of 
the  Bessel  functions  themselves  and  a method  of  hunting  for  zeros. 
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The  boundary  conditions  require 

^ = at  r = r„r, 

=*  aJm^)+bYm^)=0 
aJ'M'z)+bYm(tLtmrJ=  0 
This  only  has  non-zero  solutions  for  a,  b if 

j„,  to  K (to)=  ym  to  K to?  ) 

Using  the  recurrence  relationships 

2Jm(z)=JnM-Jm+,(z)  and  27m(z)=Tm_I(z)-ym+1(z) 


(30,31,32) 


(33) 


( 34,  35  ) 


this  equation  can  be  rewritten 

to 1 to  )-  Jm*i  »,n  )1k,-i  to?  )-  ym*i  to?  ))- 

O'— i (iJli'i )-  *'—1  (>*-'1  )X^— . )“  t-'a  ))** 0 • (36) 

Viewing  this  as  a function  of  M-l , the  values  that  are  able  to  satisfy  the  boundary  conditions  can  thus  be 

found  by  hunting  for  zeros.  The  ratio  b/a  can  be  found  by  applying  the  boundary  condition  (30)  at  either 
the  inner  or  the  outer  wall.  The  amplitude  and  phase  of  the  remaining  constant  (a  or  b ) is  arbitrary.  A 
convenient  practical  choice  is  to  set  the  value  at  the  outer  wall  to  unity,  Bkm (row ) = 1 . In  this  way,  the 

coefficient  of  the  pressure  represents  the  value  that  would  be  measured  by  a microphone  at  the  outer  wall 
in  the  absence  of  other  radial  harmonics. 


